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14.  ABSTRACT 

A  novel  wavelet  adaptive  multilevel  representation  (WAMR)  algorithm  combined  with  an 
extended  intrinsic  low  dimensional  manifold  (ILDM)  approach  to  solve  challenging 
multi scale  problems  has  been  developed.  The  emphasis  of  this  work  was  the  development 
of  a  combined  algorithm  to  rationally  reduce  to  lower  dimensions  the  differential 
equations  that  model  unsteady,  multi -dimensional  combustion  processes.  In  particular, 
there  was  a  focus  on  problems  in  which  there  is  full  coupling  between  fluid  mechanics 
and  chemical  kinetics.  Such  problems  are  characterized  by  a  multiplicity  of  phenomena 
evolving  over  length  and  time  scales  which  are  widely  disparate.  Results  are  shown  for 
a  resolved  viscous  one-dimensional  detonation  in  hydrogen/ oxygen/ argon  with  detailed 
kinetics.  The  method  has  been  extended  to  multiple  dimensions,  and  results  are  shown 
for  flow  in  a  lid-driven  cavity  at  high  Reynolds  niomber. 
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Executive  Summary 

The  scope  of  this  project  has  been  the  development  of  a  combined  algorithm  based  on  a  wavelet  adaptive 
multilevel  representation  (WAMR)  method  and  a  technique  to  rationally  reduce  to  lower  dimensions  the 
differential  equations  that  model  unsteady,  multi-dimensional  combustion  processes,  known  as  the  intrinsic 
low-dimensional  manifold,  or  ILDM. 

We  focus  on  problems  in  which  there  is  full  coupling  between  fluid  mechanics  and  chemical  kinetics.  Such 
problems  are  characterized  by  a  multiplicity  of  phenomena  evolving  over  length  and  time  scales  which  are 
widely  disparate.  In  brief,  the  wavelet  decomposition  adaptively  reduces  the  total  number  of  unknowns  in 
space,  and  the  rational  reduction,  which  applies  ILDM  theory,  subsequently  reduces  the  dimension  of  the 
local  manifold  on  which  the  chemical  reactions  evolve.  The  first  step  of  the  method  utilizes  a  multilevel 
wavelet  collocation  approximation.  The  multilevel  structure  of  the  algorithm  facilitates  the  computational 
adaptation  of  the  grid  refinement  according  to  the  temporal  evolution  of  the  solution.  One  of  the  strengths  of 
the  adaptive  wavelet  collocation  technique  is  its  ability  to  allow  high  resolution  computations  only  in  regions 
where  sharp  variations  occur.  Thus,  far  fewer  unknowns  are  needed  for  wavelets  than  for  conventional  finite- 
difference  or  finite-element  techniques.  Numerical  solutions  have  shown  the  ability  of  the  wavelet  method 
to  resolve  localized  structures,  whose  positions  and  sharpness  change  both  spatially  and  temporally.  In  the 
ILDM  phase  of  the  computation,  a  local  linear  analysis  gives  the  time  scales  associated  with  the  reactions. 
Those  processes  which  occur  on  time  scales  faster  than  a  physically  motivated  cutoff  time  scale  are  allowed 
to  relax  to  partial  equilibrium;  the  remaining  processes  are  subjected  to  time  evolution.  Consequently,  the 
number  of  differential  equations  to  be  solved  is  reduced,  and  those  that  remain  are  significantly  less  stiff. 
This  reduction  results  in  a  more  efficient  numerical  method. 

We  demonstrate  the  numerical  algorithm  by  considering  a  standard  ignition  delay  problem  for  a  mixture 
of  hydrogen/oxygen/argon  in  a  shock  tube.  An  operator  splitting  method  is  used  to  describe  the  reactions 
as  a  system  of  ordinary  differential  equations  at  each  spatial  point.  The  ILDM  method  is  used  to  eliminate 
the  stiffness  associated  with  the  chemistry  by  decoupling  processes  which  evolve  on  fast  and  slow  time 
scales.  The  fast  time  scale  processes  are  systematically  equilibrated,  thereby  reducing  the  dimension  of  the 
phase  space  required  to  describe  the  reactive  system.  The  WAMR  technique  captures  the  detailed  spatial 
structures  automatically  with  a  small  number  of  basis  functions  thereby  further  reducing  the  number  of 
variables  required  to  describe  the  system.  A  maximum  of  only  300  collocation  points  and  15  scale  levels 
yields  results  with  striking  resolution  of  fine  scale  viscous  and  induction  zones.  Additionally,  the  resolution 
of  physical  diffusion  processes  minimizes  the  effects  of  potentially  reaction-inducing  artificial  entropy  layers 
associated  with  numerical  diffusion. 

To  demonstrate  the  multidimensional  capability  of  the  WAMR  algorithm,  we  present  results  for  a  classical 
benchmark  problem:  the  two-dimensional  lid  driven  cavity.  The  problem  is  studied  at  sufficiently  large 
Reynolds  numbers  so  that  solutions  using  conventional  numerical  methods  become  difficult  to  obtain  due  to 
the  increasing  presence  of  multiscale  structures. 

The  work  is  not  incremental  in  character  and  in  fact  represents  a  substantially  new  approach  which  we 
believe  has  yielded  consequent  dividends.  While  additional  generalizations  can  and  should  be  made  to  the 
algorithm,  we  nevertheless  believe  the  research  has  set  a  new  standard  in  high  performance  computing  applied 
to  combustion.  In  addition  we  believe  the  combined  novelty  and  utility  of  the  WAMR/ILDM  method  will 
have  ramifications  extending  into  the  other  areas  of  computational  physics. 

The  following  personnel  were  involved  in  the  project: 

•  Samuel  Paolucci,  Professor, 

•  Joseph  M.  Powers,  Associate  Professor, 

•  Christopher  Bowman,  Post-Doctoral  Research  Associate, 

•  Sandeep  Singh,  Graduate  Research  Assistant, 

•  Yevgenii  Rastigejev,  Graduate  Research  Assistant, 
and  the  following  publications  resulted  from  this  work: 
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1  Introduction 


Here  we  discuss  a  novel  wavelet  adaptive  multilevel  representation  (WAMR)  algorithm  combined  with  an 
extended  intrinsic  low  dimensional  manifold  (ILDM)  approach  to  solve  challenging  multiscale  problems.  The 
emphasis  of  this  work  is  the  development  of  a  combined  algorithm  to  rationally  reduce  to  lower  dimensions 
the  differential  equations  that  model  unsteady,  multi-dimensional  combustion  processes  In  particular,  there 
is  a  focus  on  problems  in  which  there  is  full  coupling  between  fluid  mechanics  and  chemical  kinetics.  Such 
problems  are  characterized  by  a  multiplicity  of  phenomena  evolving  over  length  and  time  scales  which  are 
widely  disparate.  In  brief,  the  WAMR  algorithm  adaptively  reduces  the  total  number  of  unknowns  in 
space,  and  the  rational  reduction,  which  applies  ILDM  theory,  subsequently  reduces  the  dimension  of  the 
local  manifold  on  which  the  chemical  reactions  evolve.  The  multilevel  structure  of  the  algorithm  facilitates 
the  computational  adaptation  of  the  grid  refinement  according  to  the  temporal  evolution  of  the  solution. 
One  of  the  strengths  of  the  adaptive  wavelet  collocation  technique  is  its  ability  to  allow  high  resolution 
computations  only  in  regions  where  sharp  variations  occur.  Thus,  far  fewer  unknowns  are  needed  for  wavelets 
than  for  conventional  finite-difference  or  finite-element  techniques.  Numerical  solutions  have  shown  the 
ability  of  the  wavelet  method  to  resolve  localized  structures,  whose  positions  and  sharpness  change  both 
spatially  and  temporally.  In  the  ILDM  phase  of  the  computation,  a  local  linear  analysis  gives  the  time  scales 
associated  with  the  reactions.  Those  processes  which  occur  on  time  scales  faster  than  a  physically  motivated 
cutoff  time  scale  are  allowed  to  relax  to  partial  equilibrium;  the  remaining  processes  are  subjected  to  time 
evolution.  Consequently,  the  number  of  differential  equations  to  be  solved  is  reduced,  and  those  that  remain 
are  significantly  less  stiff.  This  reduction  results  in  a  more  efficient  numerical  method.  To  demonstrate  the 
algorithm’s  capabilities,  we  present  results  for  a  shock  tube  ignition  delay  problem  and  a  two-dimensional 
high  Reynolds  number  lid-driven  cavity  problem.  We  give  below  a  brief  discussion  of  some  of  the  challenges 
that  exist  in  simulating  realistic  combustion  problems  for  which  our  algorithm  is  well  suited. 

A  typical  combustion  event  in  devices  of  practical  importance  can  be  described  as  a  highly  coupled  process 
involving  a  large  number  of  elementary  chemical  reactions  evolving  simultaneously  within  a  complex  flow 
field.  It  is  clear  in  many  applications  that  chemical  kinetics  affects  the  motion  of  the  fluid,  and  moreover, 
the  motion  of  the  fluid  affects  the  chemical  kinetics.  One  of  the  key  problems  in  theoretical  combustion 
science  is  to  generate  reliable  solutions  to  equations  which  model  these  systems.  Often  these  models  consist 
of  equations  representing  hundreds  of  elementary  reaction  steps  coupled  with  full  compressible  Navier-Stokes 
equations.  Successful  modeling  has  the  potential  to  yield  more  precise  design  and  control  strategies  for  a 
wide  variety  of  practical  devices  such  as  propulsion  systems,  internal  combustion  engines,  incinerators,  or 
furnaces. 

In  the  past  decade,  improved  computer  hardware  has  made  practical  the  numerical  modeling  of  processes 
with  both  detailed  chemical  kinetics  and  fluid  mechanics.  Examples,  such  as  laminar  diffusion  flames,  are 
discussed  by,  among  many  others,  Smooke,  1991a.  The  kinetic  schemes  involved  are  typically  quite  large, 
involving  hundreds  of  reactions  and  dozens  of  chemical  species,  cf.  Warnatz,  1992.  As  discussed  by  Maas  and 
Pope  (MP),  1992b,  typical  reaction  kinetic  schemes  can  involve  processes  which  evolve  over  widely  disparate 
time  scales,  ranging  typically  from  10“^  s  to  10^  s.  These  time  scales  often  eclipse  those  scales  typically 
associated  with  fluid  mechanical  phenomena  which  are  governed  often  by  advection,  diffusion,  and  turbulence. 
Consequently,  fluid  mechanics  problems,  which  are  often  complicated  by  simultaneous  phenomena  occurring 
on  micro-  and  macro-scales,  are  made  even  more  daunting  by  the  addition  of  chemical  reactions  which  occur 
on  both  faster  and  slower  time  scales  than  those  of  the  fluid  mechanics.  As  a  fully  converged  numerical 
solution  to  a  set  of  equations  requires  a  full  resolution  of  all  space  and  time  scales,  completely  resolved 
solutions  of  multi-dimensional  unsteady  flow  fields  using  full  reaction  kinetics  are  essentially  impossible  with 
existing  computer  hardware  and  will  be  for  some  time  to  come. 

These  difficulties  have  motivated  a  widespread  effort  in  recent  years  to  identify  ways  to  reduce  the 
complexity  of  kinetic  models  while  simultaneously  retaining  the  flavor  of  the  full  kinetic  model,  cf.  Smooke, 
1991b;  Radhakrishnan,  1991;  Prenklach,  1991;  Prenklach,  Wang,  and  Rabinowitz,  1992;  Trevino  and  Mendez, 
1992;  Goussis  and  Lam,  1992;  Buckmaster,  1993;  Lam,  1993;  Maas  and  Pope,  1992a,  1994;  Eggles  and 
DeGoey,  1995;  Schmidt,  et  ai,  1996,  or  Vlachos,  1996.  The  emphases  in  these  studies  has  been  on  problems 
which  can  be  modeled  as  systems  of  ordinary  differential  equations,  such  as  closed  homogeneous  reactors  or 
steady  wave  propagation.  The  more  challenging  problems  in  which  systems  evolve  in  both  space  and  time 
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have  received  almost  no  scrutiny.  Such  systems  are  well-suited  for  solution  via  a  novel  wavelet  algorithm 
combined  with  an  extended  ILDM  approach.  We  have  recently  made  some  inroads  to  this  end:  Rastigejev 
and  Paolucci,  1999;  Singh,  Powers,  and  Paolucci,  1999;  Singh  and  Powers,  1999;  Powers,  Paolucci,  Rastigejev, 
and  Singh  1999;  Rastigejev,  Singh,  Bowman,  Paolucci,  and  Powers,  2000;  Singh,  Rastigejev,  Paolucci,  and 
Powers,  2001;  Rastigejev  and  Paolucci,  2002.  We  discuss  these  methods  and  the  progress  made  in  the 
following  sections. 


2  General  Problem  Description 


Here  we  describe  in  a  highly  abbreviated  form  the  set  of  unsteady,  compressible,  reactive  Navier-Stokes 
equations  in  D  (Z)  =  1,2,  or  3)  spatial  variables.  These  equations  are  a  non-linear  system  of  partial 
differential  equations  in  spatial  variables  Xi  (i  =  1, . . . ,  D)  and  time  t: 


d  ,  .  d  .  . 

St  <'’>  + 


dt 


{pUi)  +  (pUiUj  +p-Tij) 


dt 


^)) 


^'j+n,{pS,t-T.,)  +  j;) 


0, 

0, 

0 

0, 


(1) 

(2) 

(3) 

(4) 

(5) 


Dependent  variables  are  density  p\  velocity  vector  Ui\  internal  energy  e;  element  mass  fractions  yi^  where 
I  =  1, . . . ,  L  —  1,  and  there  are  L  elements;  species  mass  fractions  Yn^  where  n  =  1, . , . ,  ZZ  —  L,  and  there 
are  N  species,  pressure  p,  viscous  stress  Tij,  energy  flux  J?,  element  flux  j/,  diffusive  species  flux  J^,  and 
species  production  rate  Here  M.n  is  the  species  molecular  mass,  and  Sij  is  the  Kronecker  delta.  The 
above  equations  describe  the  conservation  of  mass,  linear  momenta,  energy,  and  element  mass  fractions, 
respectively,  along  with  the  time  evolution  of  species  mass  fractions.  The  set  of  equations  is  completed  when 
constitutive  equations  for  energy,  pressure,  heat  flux,  viscous  stress,  and  species  production  rates  are  given. 
Examples  of  complete  formulation  of  these  terms  is  standard  and  given,  for  instance  by  Buckmaster  and 
Ludford,  1982,  or  Williams,  1985. 

The  governing  equations  can  be  written  in  the  form 


= -■^f{u{xi,t))  +  g{u{xi,t)),  u,f,g  €  (6) 

u=  (p,pui,p(e  +  ^^^^  ,pyi,pYn'j  ■  (7) 

Here,  we  have  defined  u  as  a  vector  of  length  -f  D  -f  1  dependent  variables.  The  vector  f,  which  can 
be  written  as  a  function  of  u,  models  the  advective  and  diffusive  fluxes  of  u,  and  the  vector  g  represents  a 
source  due  to  chemical  reaction. 


3  Operator  Splitting 


Equation  (6)  is  solved  in  two  steps,  a  reaction  step  and  a  convection-diffusion  step,  using  Strang-splitting 
(Strang,  1968)  .  This  splitting  results  in  second  order  accuracy  in  time.  The  set  of  partial  differential 
equations  is  reduced  to  ordinary  differential  equations  at  each  spatial  node  /c  by  a  wavelet  collocation  method 
to  be  described  shortly.  In  the  first  step,  we  suppress  advection  and  diffusion,  and  account  for  reaction. 
Suppression  of  advection  and  diffusion  results  in  homogeneous  equations  for  mass,  linear  momenta,  energy, 
and  element  maiss  fraction  which  can  be  integrated  analytically  so  as  to  leave  only  N  —  L  equations  to 
integrate.  We  then  suppress  reaction  and  solve  the  ODEs  which  account  for  advection  and  diffusion  using 
the  wavelet  method. 


4 


3.1  Reaction  step 


In  Step  1,  the  following  equation,  equivalent  to  that  for  a  homogeneous  premixed  reactor  at  every  point 
in  the  space,  is  solved  at  each  spatial  point 

|=g(u).  (8) 

The  first  1  +  D  +  L  equations  of  Equation  (8)  are  homogeneous  and  can  be  integrated  exactly  to  give 


p  =  po  Ui^Uio  e=^eo  yi^yio  (/ =  1). 


(9) 


Hence,  the  N  —  L  species  evolution  equations  in  Equation  (8)  reduce  to  the  following  partial  differential 
equations,  which  are  treated  at  each  point  in  space  as  ordinary  differential  equations  with  p,  Wi,  e,  and  yi 
held  constant  to  the  values  given  in  Equation  (9)  every  time  this  step  is  repeated: 


OYj  ^  CjjMj 
dt  po 


(10) 


It  is  noted  that  during  a  single  time  step  the  values  of  Po^^o^  and  yio  do  not  change;  because  of  convection 
and  diffusion,  in  general,  all  will  change  with  time. 

It  can  be  shown  that  Ui  can  be  cast  in  terms  of  for  z  —  1, . . . ,  iV  — L;  thus,  Equations  (10)  are  well  posed 
and  can  be  solved  in  their  entirety  by  any  standard  implicit  or  explicit  technique  at  every  point  in  space. 
Because  the  equations  are  stiff,  we  can  use  the  LSODE  software  package,  Hindmarsh,  1983,  in  full  implicit 
mode  when  solving  the  full  set  of  Equations  (10).  Alternatively,  Equations  (10)  can  be  solved  using  the  ILDM 
method,  described  in  a  following  section,  which  systematically  removes  the  stiffness  associated  with  reactions. 
When  using  the  ILDM  method,  we  use  a  simpler  explicit  Runge-Kutta  method  for  time  advancement,  which 
is  second  order  accurate  in  time.  The  size  of  the  time  step  is  dictated  by  the  convection-diffusion  time  step 
restriction. 


3.2  Convection-diffusion  step 


In  Step  2,  the  following  equation,  which  is  a  set  of  partial  differential  equations,  is  solved  for  the  convection 
diffusion  step: 


0u 

dt 


(11) 


Again,  Equation  (11)  can  be  solved  by  a  variety  of  standard  discretization  techniques  developed  for  inert 
fluid  mechanics.  Here  we  use  the  WAMR  technique,  as  discussed  in  detail  by  Vasilyev  and  Paolucci,  1996, 
1997.  The  method  is  summarized  as  follows.  At  any  given  time  step,  the  pressure,  temperature,  density,  and 
velocity  fields  are  projected  onto  a  multilevel  wavelet  basis.  The  amplitudes  of  the  wavelet  basis  functions 
give  a  measure  of  the  importance  of  a  particular  wavelet  mode.  Additionally,  one  has  available  a  priori 
error  estimates,  in  contrast  to  most  gradient-based  adaptive  mesh  refinement  techniques.  All  wavelets 
whose  amplitude  are  below  a  defined  threshold  are  removed.  Calculations  are  performed  for  each  wavelet 
whose  amplitude  is  above  the  threshold  (essential  wavelets)  and  for  a  certain  number  below  the  threshold 
(neighboring  wavelets).  If  at  the  completion  of  a  time  step,  an  essential  wavelet  has  its  amplitude  drop 
below  the  threshold,  it  is  reclassified  as  a  neighboring  wavelet  or  eliminated,  and  the  neighboring  region 
is  adjusted;  similarly,  if  a  wavelet  in  the  neighboring  region  has  its  amplitude  become  sufficiently  large,  it 
is  reclassified  as  an  essential  wavelet,  and  the  neighboring  region  is  adjusted.  The  method  is  based  on  a 
collocation  strategy  using  the  auto- correlation  functions  of  the  Daubechies  scaling  function  and  wavelet  of 
order  four  as  bases.  A  linearized  trapezoidal  (implicit)  scheme  in  conjunction  with  GMRES  iterations  is 
used  for  time  advancement.  The  size  of  the  time  step  is  chosen  to  satisfy  a  CFL  condition  associated  with 
the  fastest  local  velocity. 
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4  Intrinsic  Low  Dimensional  Manifolds 

As  noted  above,  each  discrete  point  in  space  behaves  as  an  adiabatic,  homogeneous  premixed  reactor, 
thus  enabling  the  use  of  the  ILDM  technique.  The  method  provides  a  systematic  way  to  overcome  the 
severe  stiffness  which  is  associated  with  full  models  of  gas  phase  combustion,  and  thus  significantly  improves 
computational  efficiency.  The  ILDM  method  systematically  eliminates  the  stiffness  by  equilibrating  fast  time 
scale  events  and  describing  parametrically  a  low  dimensional  manifold  upon  which  slow  time  scale  events 
evolve. 

Here,  the  method  of  Maaa  and  Pope,  1992b,  slightly  adapted  for  our  system,  is  summarized.  The 
constant  density,  isochoric,  adiabatic  combustion  for  the  homogeneous  premixed  system  is  described  by 
Eq.  (10).  The  remaining  L  mass  fractions  are  computed  from  algebraic  constraints  described  in  full  by 
Singh,  et  a/.,  2001.  For  gases  in  the  mixture,  which  are  thermally  perfect,  the  specific  internal  energy  of  each 
component  is  at  most  a  function  of  temperature.  Hence,  knowledge  of  the  mass  fractions  and  mass-averaged 
specific  internal  energy  Cq  allows  one  to  use  Newton’s  method  to  invert  the  caloric  state  equation  to  form 
T  =  T(yi, . . . ,  IW-l;  Co,  yic )  •  •  • » y(L-i)o)*  Hence,  which  is  in  general  a  function  of  temperature,  density, 
and  species  mass  fractions,  can  now  be  considered  as  a  function  of  only  —  L  of  the  species  mass  fractions 
for  fixed  values  of  e^,  /?o,  and  yio'. 

Ui  ==a;i(yi,...,y}v-L;eo,/9o,yu,...,y(L-i)o)>  ^  =  l,...,A/’-L.  (12) 

Both  uJi  and  et  can  be  easily  evaluated  using  the  Chemkin  HI,  Kee,  et  a/.,  1990,  package.  Equation  (10)  can 
thus  be  rewritten  as 

^  -  G(Y),  (13) 

where  Y  =  (Yu...,Yn-lV  and  G(Y)  =  . . . , The  reciprocals  of  the  eigenvalues  of 

the  Jacobian  Gy  at  any  point  in  the  composition  space  identify  the  N  —  L  associated  characteristic  time 
scales.  The  corresponding  eigenvectors  identify  the  directions  in  which  each  eigenmode  of  the  total  trajectory 
evolves  in  composition  space.  The  ILDM  is  identified  by  a  set  of  points  in  the  composition  space  where  the 
composition  space  velocity  vector  G  is  orthogonal  to  the  eigenvectors  associated  with  fast  time  scales. 

We  outline  here  the  analysis  to  compute  the  ILDM  given  in  detail  for  a  similar  system  first  by  Maas 
and  Pope,  1992b,  and  later  extended  by  Maas,  1998,  following  the  general  theory  given  by  Golub  and  Van 
Loan,  1996.  We  choose  to  resolve  M  slow  time  scales;  this  will  result  in  the  choice  of  M  species  which  will 
form  part  of  the  parameterization  for  the  complete  ILDM.  This  is  formed  via  the  following  steps.  First  we 
decompose  the  Jacobian  matrix  Gy  with  a  real  Schur  decomposition: 

Gy  =  Q-U.Q^,  (14) 

where  Q  is  the  {N  —  L)  x  {N  —  L)  orthogonal  matrix  of  real  Schur  vectors  and  U  is  an  {N  —  L)  x  {N  —  L) 
upper-triangular  matrix  with  eigenvalues  along  the  main  diagonal.  Typically  most  eigenvalues  are  negative 
and  real;  for  closed  adiabatic  systems  near  equilibrium,  there  is  a  guarantee  of  real  negative  eigenvalues.  The 
matrix  is  ordered  so  that  proceeding  down  the  diagonal,  the  real  part  of  the  eigenvalue  becomes  progressively 
more  negative.  Thus,  the  eigenvalue  at  the  bottom  right  is  associated  with  the  fastest  stable  mode. 

The  above  Schur  decomposition,  combined  with  the  solution  of  a  Sylvester  equation  allows  the  Jacobian 
to  be  written  in  the  following  form 

Gv  =  (Zs  N,)(lp) 

Here  Z5  has  M  column  vectors,  each  of  length  N  -  L,  which  span  the  subspace  associated  with  the  M  slow 
time  scales.  Further,  Zf  has  V  -  L  -  M  column  vectors,  each  of  len^h  AT  ~  L,  which  span  the  subspace 
associated  with  the  N  —  L~  M  fast  time  scales.  The  matrices  Z5  and  Zf,  as  defined  in  Equation  (15),  have 
dimension  Af  x  (iV  —  L)  and  (iV  —  jL  —  A/)  x  (A/^  —  L),  respectively,  and  are  associated  with  the  reciprocal  slow 
and  fast  basis  vectors.  The  matrix  N5  has  dimension  M  x  M,  is  upper- triangular,  and  has  as  its  eigenvalues 
those  associated  with  the  slow  time  scales.  The  matrix  Nf  has  dimension  {N  ~  L  —  M)  x  {N  —  L  —  M),  is 
upper-triangular,  and  has  as  its  eigenvalues  those  associated  with  the  fast  time  scales. 
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For  fixed  energy,  density,  and  element  mass  fractions,  the  ILDM  method  identifies  M-dimensional  sub¬ 
spaces  (M  <  N  —  L)  on  which  slower  time  scale  events  evolve.  An  ILDM  which  equilibrates  the  N  —  L  —  M 
fast  time  scales  and  consequently  describes  the  M-dimensional  manifold  is  given  by 

Zf  •  G(Y)  -  0,  (16) 

which  forms  N  —  L  —  M  algebraic  equations.  For  what  is  effectively  a  complicated  one-step  chemistry,  we 
take  M 

To  construct  the  manifold  for  fixed  energy,  density,  and  element  mass  fractions,  one  first  determines  the 
equilibrium  point  of  the  system,  which  is  taken  as  the  original  seed  value.  In  subsequent  calculations,  the 
most  recently  calculated  point  is  used  as  the  seed  value.  One  then  performs  the  local  fast  and  slow  subspace 
decomposition,  which  gives  local  eigenvalues  and  associated  basis  vectors  for  a  system  linearized  about  the 
seed  value.  One  then  perturbs  —  L  of  the  species  away  from  their  seed  values  in  a  prescribed  manner  to 
form  the  M  algebraic  equations  below: 

P.(Y-Yo)  =  <5.  (17) 

Here  P  is  a  user-specified  parameterization  matrix  of  dimension  M  x  (A^  —  L),  5  is  a  vector  of  length  M 
which  contains  a  user-specified  increment  in  the  projected  mass  fractions,  and  Yo  is  a  vector  of  length 
N  —  L  which  contains  the  seed  values  of  mass  fractions.  These  M  equations  are  solved  simultaneously  with 
the  N  —  L  —  M  algebraic  ILDM  equations  (16)  with  a  predictor  corrector  technique  which  uses  a  tangent 
predictor  and  a  Newton’s  method  corrector  to  obtain  the  mass  fractions  at  a  new  point  on  the  manifold. 
The  process  is  repeated  to  construct  the  global  ILDM.  We  adapt  the  parameterization  matrix  P  to  overcome 
difficulties  associated  with  turning  points  of  the  manifold  in  composition  space.  The  choice  of  P  is  based 
on  a  continuation  method  which  is  an  extrapolation  of  previous  points  of  the  ILDM  used  to  estimate  the 
location  of  the  new  points.  This  method  can  also  be  used  for  computing  manifolds  for  an  adiabatic,  isobaric 
system. 

With  this  analysis,  an  M-dimensional  manifold  can  be  identified  in  an  {N  ~  L)-dimensional  composition 
space  for  a  given  set  of  po?  ^  Thus,  a  different  ILDM  is  required  for  a  different  set  of 

densities,  internal  energies  and  element  mass  fractions.  Since  in  general  calculations  one  can  expect  all 
of  these  quantities  to  vary,  the  actual  relevant  reaction  manifold  which  must  be  formed  has  dimension 
K  =  M  1  L  and  can  be  tabulated  numerically  to  give 

yi  =  Yi{Yu...,YM,e,p,yi,...,yL-i)  for  i  =  M  +  I, . . .  ,N  -  L,  (18) 

where  Yi, . . . ,  Ym  are  the  chosen  reference  variables  for  the  ILDM  lookup  table.  The  reference  variables  are 
chosen  in  such  a  way  that  the  ILDM  is  single  valued  with  respect  to  these  variables  for  easy  lookup.  While 
there  is  no  guarantee  of  single- valuedness,  in  the  problems  studied  we  have  found  it  to  be  the  case. 

It  is  the  dimension  of  the  lookup  table,  AT,  which  is  critically  important  in  manifold  methods.  Many 
previously  reported  calculations  have  been  restricted  to  premixed  conditions  in  the  isobaric  and/or  adiabatic 
limits,  thus  reducing  the  dimension  K.  In  the  calculations  presented  here,  variable  energy  and  density  are 
intrinsic  features  of  the  flow;  given  that  we  have  chosen  M  =  1  for  the  hydrogen /oxygen/ argon  system 
discussed  later  (L  =  3)  we  in  principle  must  deal  with  a  table  which  has  dimension  K  ~  5.  Because  we 
will  study  a  uniformly  premixed  mixture,  and  because  we  employ  a  simplified  diffusion  model  in  which 
there  is  no  preferential  diffusion,  L  —  1  =  2  independent  element  mass  fractions  remain  constant  throughout 
the  calculations,  reducing  the  effective  dimension  of  the  manifold  to  AT  =  3.  The  variables  of  manifold 
parameterization  for  our  hydrogen/argon/oxygen  problem  are  chosen  to  be  e,  p,  and  Y//20-  In  the  ignition 
delay  problem  considered  later  for  a  mixture  of  hydrogen/  oxygen/argon,  we  use  a  9  species  37  step  reaction 
mechanism  described  in  Table  1.  With  M  =  1,  a  projection  of  the  ILDM  for  fixed  internal  energy,  density, 
(as  well  as  the  always  fixed  element  mass  fractions),  for  this  system  is  plotted  in  Fig.  1  with  Y//2O  nsed  as 
the  reference  independent  variable  for  the  ILDM  and  YH2O2  ^  the  dependent  variable.  Mass  fractions  for 
all  species  (not  shown  here)  are  also  available  as  functions  of  Y//2O  for  the  same  fixed  internal  energy  and 
density.  Also  shown  on  the  plot  are  projections  of  trajectories  in  this  two-dimensional  subspace  for  a  variety 
of  initial  conditions.  It  is  seen  that  all  trajectories  relax  to  the  ILDM.  As  indicated  on  Fig.  1,  the  relaxation 
from  the  initial  state  to  the  manifold  occurs  on  a  relatively  fast  time  scale  on  the  order  of  microseconds, 
while  once  on  the  ILDM,  the  subsequent  relaxation  to  final  equilibrium  occurs  on  a  much  slower  time  scale 
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Reaction 

aj 

/Si 

E, 

1 

02  +  H^0H  +  0 

2.00  X  10^"* 

0.00 

70.30 

2 

OH  +  O  -^02  +  H 

1.46  X  IQi^ 

0.00 

2.08 

3 

H2  0  — >  OH  H 

5.06  X  lO" 

2.67 

26.30 

4 

OH  +  H  H2  +  O 

2.24  X  10^ 

2.67 

18,40 

5 

H2  +  OH  ^  H2O  +  H 

1.00  X  108 

1.60 

13.80 

6 

H2O  +  H-^H2  +  0H 

4.45  X  10® 

1.60 

77.13 

7 

OH  +  OH  ->  H2O  +  0 

1.50  X  10® 

1.14 

0.42 

8 

H2O  +  0-*0H  +  0H 

1.51  X  lOi® 

1.14 

71.64 

9 

H  +  H  +  M-*H2  +  M 

1.80  X  10^® 

-1.00 

0.00 

10 

H2  +  M^H  +  H  +  M 

6.99  X  10^® 

-1.00 

436.08 

11 

H  +  OH +  M  ^  H2O  +  M 

2.20  X  10^2 

-2.00 

0.00 

12 

H2O  +  M^H  +  OH  +  M 

3.80  X  1023 

-2.00 

499.41 

13 

0  -f-  0  “h  M  — >  O2  4“  M 

2.90  X  10*’’ 

-1.00 

0.00 

14 

O2  “h  A/  — >  0  “h  0  “h  Af 

6.81  X  10‘® 

-1.00 

496.41 

15 

H  4“  O2  4“  A/  — ^  HO2  4~  Af 

2.30  X  lOi® 

^0.80 

0.00 

16 

HO2  +  M  H  +  O2  +  M 

3.26  X  lOi® 

-0.80 

195.88 

17 

HO2  +  H  ^  OH +  OH 

1.50  X  10‘^ 

0.00 

4,20 

18 

OH  +  OH  HO2  +  H 

1.33  X  1013 

0.00 

168.30 

19 

HO2  +  H-*H2  +  02 

2.50  X  1013 

0.00 

2.90 

20 

H2  +  02-^  HO2  +  H 

6.84  X  1013 

0.00 

243.10 

21 

HO2  +  H  H2O  +  0 

3.00  X  1013 

0.00 

7.20 

22 

H2O  +  0  HO2  +  H 

2.67  X  1013 

0.00 

242.52 

23 

H  O2  +  0  — *  OH  +  O2 

1.80  X  1013 

0.00 

-1.70 

24 

OH +  02^  HO2  +  0 

2.18  X  1013 

0.00 

230.61 

25 

HO2  +  OH  -*  H2O  +  O2 

6.00  X  1013 

0.00 

0.00 

26 

H2O  +  O2-*  HO2  +  OH 

7.31  X  1014 

0.00 

303.53 

27 

HO2  +  HO2  H2O2  +  O2 

2.50  X  lOii 

0.00 

-5.20 

28 

OH  +  OH  +  M  ^  H2O2  +  M 

3.25  X  1022 

-2.00 

0.00 

29 

H2O2  +  M  OH  +  OH  +  M 

2.10  X  1024 

-2.00 

206.80 

30 

H2O2  +  H^H2  +  HO2 

1.70  X  IO12 

0.00 

15.70 

31 

H2  +  HO2  ^  H2O2  +  H 

1.15  X  IO12 

0.00 

80.88 

32 

H2O2  +  H^  H2O  +  OH 

1.00  X  1013 

0.00 

15.00 

33 

H2O  +  OH  ^  H2O2  +  H 

2.67  X  IO12 

0.00 

307.51 

34 

H2O2  +  0  ^  OH  +  HO2 

2.80  X  1013 

0.00 

26.80 

35 

OH  +  HO2  ^  H2O2  +  0 

8.40  X  IO12 

0.00 

84.09 

36 

H2O2  +  OH  ^  H2O  +  HO2 

5.40  X  IO12 

0.00 

4.20 

37 

H2O  +  HO2  H2O2  +  OH 

1.63  X  1013 

0.00 

132.71 

Table  1:  Nine  species,  thirty-seven  step  reaction  mechanism  for  hydrogen/oxygen/argon  mixture,  Maas  and 
Warnatz,  1988,  with  corrected  ///2  from  Maas  and  Pope,  1992,  also  utilized  by  Fedkiw,  et  a/.,  1997.  Units 
of  aj  are  in  appropriate  combinations  of  cm,  moZ,  s,  and  K  so  that  Coi  has  units  of  mol  cm  ^  s  units  of 
Ej  are  kJ  mol~^.  Third  body  collision  efficiencies  with  M  are  ///2  =  1.00,  fo2  =  0.35,  and  ///2O  =  6.5. 
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Figure  1:  ILDM  projection  for  9  species  37  step  reaction  mechanism  of  hydrogen/oxygen/argon  combustion 
as  function  of  YH2O  at  constant  values  of  p  =  0.5  kg  e  =  8.0  x  10^  J  Element  mass  fractions 

fixed  at  yn  =  0.01277,  yo  =  0.10137,  =  0.88586.  Also  plotted  are  trajectories  from  full  time  integration 

showing  relaxation  to  the  manifold  and  equilibrium.  The  symbol  ‘x’  denotes  equally  spaced  0.10  ps  time 
intervals.  Total  time  to  relax  to  equilibrium  is  near  0.1  ms. 


on  the  order  of  0.1  ms.  The  phase  trajectories  appear  to  cross,  but  this  is  because  they  have  been  projected 
into  a  lower  dimensional  space. 

A  projection  of  the  K  ~  3  ILDM  for  constant  density  and  the  same  constant  element  mass  fractions  for 
the  hydrogen/oxygen/argon  system  is  plotted  in  Fig.  2,  with  1^2  0  and  e  used  as  the  reference  independent 
variables  and  Yh^o2  as  the  dependent  variable.  The  portion  of  the  ILDM  depicted  closely  corresponds 
to  values  realized  in  the  ignition  delay  problem  discussed  later.  Also,  for  the  system  studied  here,  the 
dependency  of  all  variables  on  density  was  weak  for  the  range  of  realized  density  values. 

Equation  (13)  can  be  solved  using  the  ILDM  method  if  the  composition  of  the  mixture  at  a  given  spatial 
point  is  close  to  the  ILDM.  One  illustrative  way  to  achieve  this  is  to  locally  project  Equation  (13)  onto  the 
slow  subspace: 

Zs-jY  =  Zs-  G(Y).  (19) 

This  results  in  M  ordinary  differential  equations  and  models  the  dynamics  of  reduced  chemistry.  In  practice, 
we  simply  integrate  M  of  the  species  equations  in  their  original  form  while  using  the  manifold  to  obtain  all 
necessary  secondary  variables. 

When  using  ILDMs  with  M  =  1,  we  need  to  integrate  only  one  differential  equation  associated  with 
the  reference  variable  at  each  spatial  location.  If  the  composition  is  far  from  the  ILDM,  it  is  important  to 
incorporate  the  off-manifold  kinetics  in  some  fashion.  Here  the  full  reaction  kinetic  equations  are  integrated 
using  LSODE  in  implicit  mode.  We  switch  from  full  integration  to  the  ILDM  method  when  the  L2  norm 
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Figure  2:  ILDM  projection  for  9  species  37  step  mechanism  of  hydrogen/oxygen/argon  combustion  giving 
FH2O2  ^  ^  function  oiYji^o  and  e  for  constant  p  —  0,5  kg  Element  mass  fractions  fixed  at  yn  =  0.01277, 
yo  =  0.10137,  yAr  =  0.88586. 


between  the  actual  mass  fractions  at  a  particular  location  and  the  corresponding  projected  mass  fractions  is 
less  than  1  x  10"®.  While  this  reduces  the  efficacy  of  the  ILDM  method,  it  is  critical  to  avoid  large  phase 
errors  associated  with  projecting  onto  the  manifold  from  a  remote  region  of  phase  space.  In  other  words, 
while  all  processes  are  typically  destined  to  reach  the  manifold,  it  is  critical  for  the  proper  sequencing  of 
events  that  they  reach  the  manifold  at  the  correct  time,  and  reach  the  right  point  on  the  manifold.  A  naive 
projection  can  easily  give  plausible  answers  which  have  0(1)  errors  associated  with  them.  In  practice,  in  our 
calculations  for  the  ignition  delay  problem,  we  find  we  are  able  to  use  the  ILDM  method  in  cells  which  have 
recently  been  shocked,  have  passed  through  the  bulk  of  the  induction  zone,  and  are  well  within  the  thermal 
explosion  region.  The  ILDM  method  is  also  applied  to  points  in  the  trailing  rarefaction  wave. 

The  convection-diffusion  step  is  equivalent  to  a  perturbation  off  the  ILDM.  Subsequent  to  the  perturba¬ 
tion,  there  is  a  fast  relaxation  to  the  manifold  corresponding  to  a  new  set  of  conserved  variables.  This  is 
accomplished  here  by  projection  onto  the  manifold  while  holding  the  reference  variable,  here  YfraOi  constant. 
This  projection  is  allowed  by  the  large  time  scale  difference  between  slow  chemical/fiuid  time  scales  and 
fast  chemical  time  scales  as  long  as  convective  and  diffusive  effects  are  not  large.  We  also  note  that  for  this 
hydrogen /oxygen/argon  mechanism  the  fast  directions  are  nearly  orthogonal  to  the  slow  variable,  YH2O} 
the  ILDM.  A  slightly  more  accurate  alternative  would  be  to  project  to  the  manifold  in  the  direction  of  the 
vectors  associated  with  the  fast  eigenmodes.  In  contrast  to  many  implementations  of  the  ILDM  method 
for  partial  differential  equations,  which  are  often  confined  to  a  single  adiabatic,  isobar ic  ILDM,  we  convect 
and  diffuse  all  variables,  not  just  slow  variables  associated  with  the  ILDM.  While  this  comes  at  a  cost  of 
solving  more  equations,  it  is  necessary  to  preserve  the  consistency  of  the  Strang  operator  splitting  method. 
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Alternatively,  for  our  problem,  which  is  neither  adiabatic,  isobaric,  or  isochoric,  we  could  choose  to  incur 
the  extra  expense  of  tabulating  the  slow  and  fast  subspaces  associated  with  changes  in  energy  and  density 
and  then  convect  and  diffuse  just  slow  variables. 

5  Wavelet  Adaptive  Multilevel  Representation  (WAMR) 


Much  of  the  material  of  this  section  builds  on  Vasilyev  and  Paolucci,  1995,  1996,  1997  and  is  described 
more  fully  in  Rastigejev  and  Paolucci,  1999,  2002. 

5.1  Bases  and  Wavelet  Approximation 

Wavelet  analysis  is  a  new  numerical  concept  (c/.  Daubechies,  1992;  Meyer,  1990),  which  allows  one  to 
represent  a  function  in  terms  of  basis  functions,  called  wavelets,  which  are  localized  in  both  location  and  scale. 
The  Daubechies  autocorrelation  functions  and  boundary-modified  Daubechies  autocorrelation  functions  (see 
Saito  and  Beylkin,  1993)  of  order  p  have  been  used  as  elements  for  the  construction  of  basis  functions.  If 
(p{x)  is  the  Daubechies  scaling  function  of  order  p  {p  coefficients  in  the  dilatation  equation  and  corresponding 
to  p/2  vanishing  moments)  then  (j){x)  =  Is  the  Daubechies  autocorrelation  function  of 

order  p  having  p  vanishing  moments.  In  addition,  following  Donoho,  1992,  we  use  the  following  definition 
for  the  wavelet: 

i){x)  =  <t>{2x-l),  (20) 

For  this  wavelet  basis,  the  algorithm  for  finding  wavelet  expansion  coefficients  reduces  to  the  use  of  simple 
Lagrange  interpolation  (see  Saito  and  Beylkin,  1993,  and  Dubuc,  1986). 

Now  any  continuous  function  f{x)  defined  in  a  finite  domain  may  be  approximated  within  a  prescribed 
accuracy  (for  L  large  enough)  by 


L 

P^/(x)  =Y^fk<j>ix-k)  +  Y2'^di,ki’i‘2‘x-k).  (21) 

k  1=0  k 

Left  and  right  boundary-modified  scaling  functions  can  be  constructed  to  span  the  finite  domain  but  we  do 
not  write  them  here  (see  Rastigejev  and  Paolucci,  1999). 

5.2  Sparse  Wavelet  Approximation  and  Irregular  Grid 


The  spatially  adaptive  collocation  grid  is  especially  well  tailored  for  strongly  time-dependent  problems. 
The  method  incorporates  the  dynamically  adaptive  multilevel  algorithm  of  Liandrat  and  Tchamitchian, 
1990.  The  adaptation  is  achieved  by  retaining  only  those  wavelets  whose  amplitudes  are  greater  than  a  given 
threshold.  Thus,  high  resolution  computations  are  carried  out  only  in  those  regions  where  sharp  transitions 
occur.  In  order  to  make  the  wavelet  approximation  of  a  function  efficient  (z.e.  to  use  the  smallest  possible 
number  of  wavelet  modes)  while  at  the  same  time  retaining  the  accuracy  of  the  representation,  we  introduce 
the  sparse  wavelet  representation  by  thresholding  the  wavelet  coefficients.  The  approximation  can  be  split 
into  two  parts  which  are  composed  of  wavelets  whose  amplitudes  are  above  and  below  a  specified  threshold 


where 


P^f{x)  =  P^f{x)  +  P^f{x), 

L 

k  i=0  fc,|d,,k|>€ 

L 

p</w  =  I]  Y.  di,kipi,k{^)- 

1=0  k,\di^k 


(22) 

(23) 

(24) 
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We  now  call  P^f{x)  the  sparse  wavelet  representation  of  the  function  /(x),  and  we  call  the  corresponding 
wavelets  used  in  the  representation  essential  wavelets.  In  this  representation,  if  a  wavelet  is  omitted,  so  is 
the  corresponding  collocation  point.  This  compression  results  in  the  irregular  grid  G>. 

5.3  Time  Integration  Algorithm 

After  the  operator  splitting  and  application  of  the  wavelet  discretization  in  space,  Equation  (11)  reduces 
to  a  system  of  ordinary  differential  equations  (ODEs): 

f = 

where  u  is  a  vector  of  dependent  variables  defined  on  the  nonuniform  collocation  grid,  and  F  represents 
the  discretization  of  the  convection  and  diffusion  terms.  Spatial  derivatives  are  approximated  by  use  of  the 
Lagrange  polynomial  corresponding  to  the  wavelet  basis  evaluated  at  collocation  points. 

To  integrate  this  system,  we  use  the  linearized  trapezoidal  method  to  obtain  the  following  system  of 
linear  algebraic  equations: 

(26) 

where  A(”)  =  (I  -  ^  [F  +  F  I  is  the 

identity  matrix  and  J  =  dF/dn  is  the  Jacobian  matrix. 

In  order  to  take  advantage  of  the  sparseness  of  matrix  as  well  as  the  fact  that  is  small,  the 

system  of  algebraic  equations  is  solved  by  using  a  multiscale  Schwarz- type  domain-decomposition  method 
which  utilizes  a  multigrid-type  iteration  cycle  for  solving  the  above  linear  problem.  First  the  hierarchy  of 
grids  are  projected  to  the  physical  domain  and  they  are  subsequently  decomposed  into  multiple  disjoint 
subsets  consisting  of  grids  of  different  size  as  shown  in  Figure  3.  The  next  step  consists  of  the  solution  of 

G,  G,  G,  G3  G. 

(• - }-<• - .J - h - - * 


Figure  3:  Iteration  procedure. 


the  linear  system  on  these  disjoint  domains  using  a  multigrid-type  iteration  cycle  on  every  subset  (linear 
subsystem).  Since  the  problem  is  often  stiff  on  the  finer  grids,  the  GMRES  method  is  used  as  the  basic 
iteration  algorithm  (see  Saad  and  Schultz,  1986).  The  method  has  potential  for  additional  development  and 
improvement. 

A  typical  cycle  is  illustrated  in  Figure  4  where  Gi  and  G4  correspond  to  the  finest  and  coarsest  grid 
subsets,  respectively.  The  number  of  unknowns  in  each  subset  is  approximately  the  same,  so  that  in  this 
example  each  subset  would  contain  approximately  one-fourth  of  the  total  number  of  grids.  Typically  the 
cycle  starts  from  the  finest  grid  subset.  The  residual  tolerance  is  decreased  on  every  successive  branch  of  the 
cycle.  We  note  that  the  number  of  iterations  required  to  satisfy  the  tolerance  drops  rapidly.  To  improve  the 
performance,  a  reasonable  (but  not  very  accurate)  initial  guess  of  the  solution  is  used.  Two  alternatives  have 
been  tried.  In  one,  the  initial  guess  is  obtained  by  solving  exactly  the  linear  tridiagonal  system  resulting 
from  a  three  point  finite-difference  derivative  spatial  approximation.  In  the  other,  we  use  the  initial  guess 


Au(^+^)  = 


{ 


Au(^), 

/(Au(^+i)), 


if  the  grid  has  not  changed, 
if  the  grid  has  changed, 


(27) 


where  we  define  Au^®^  =  0,  and  I  is  the  wavelet  interpolation  operator  which  interpolates  from  the  nonuni¬ 
form  old  grid  to  the  nonuniform  new  grid. 


5.4  Dynamically  Adaptive  Algorithm 
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Figure  4:  Iteration  procedure. 


In  order  to  resolve  all  structures  appearing  in  the  solution  of  the  problem  as  time  evolves,  we  need 
to  adapt  the  sparse  wavelet  representation  in  time.  To  accomplish  this,  at  any  time  we  need  to  save  the 
amplitudes  of  essential  wavelets  and  the  amplitudes  of  wavelets  in  the  neighboring  region^  i.  e.  those  wavelets 
which  are  near  the  essential  wavelets  in  location  and  scale.  Essential  and  neighboring  wavelets  constitute 
the  group  of  active  wavelets.  The  irregular  grid  G>  is  composed  of  collocation  points  corresponding  to  all 
active  wavelets  at  any  given  moment  of  time. 

The  numerical  algorithm  for  solving  the  system  of  ODEs  may  be  summarized  as  follows: 

1.  integrate  the  ODEs  to  obtain  the  approximate  solution  on  the  irregular  grid  by  using 

as  initial  conditions; 

2.  obtain  the  new  grid  based  on  the  magnitudes  of  wavelet  coefficients  of  and  accounting 

for  the  new  neighboring  region; 

3.  if  =  G^^  we  increment  time  and  go  back  to  step  1,  otherwise  compute  the  values  of  at 

the  new  collocation  points  in  then  increment  time  and  go  back  to  step  1. 

The  dynamically  adaptive  wavelet  collocation  algorithm  is  ideally  suitable  for  handling  problems  with 
localized  structures,  which  might  occur  intermittently  anywhere  in  the  computational  domain  or  change 
their  locations  and  scales  in  space  and  time.  The  implementation  of  conventional  adaptive  algorithms  is  very 
costly  due  to  the  possibility  of  grids  varying  drastically  within  a  short  time  interval.  The  use  of  conventional 
algorithms  on  a  uniform  grid  is  impractical.  Thus,  the  main  advantage  of  the  dynamically  adaptive  wavelet 
collocation  algorithm  is  that  it  requires  a  far  fewer  number  of  unknowns  than  the  other  algorithms,  when 
applied  to  problems  with  a  great  range  of  spatial  scales.  In  addition,  the  computational  grid  can  be  refined 
locally  to  an  arbitrary  small  size.  We  emphasize  here  that  the  adaptation  of  the  computational  grid  does 
not  require  additional  efforts  and  consists  merely  in  turning  on  and  off  wavelets  at  different  locations  and 
scales.  Moreover,  the  algorithm  exhibits  an  exponential  convergence  on  the  unevenly  spaced  mesh.  The 
computational  cost  is  independent  of  the  spatial  dimensionality  and  is  G(7V),  where  jVis  the  total  number 
of  collocation  points.  This  computational  efficiency  is  possible  by  virtue  of  the  Fast  Wavelet  Collocation 
Transform  (see  Vasilyev  and  Paolucci,  1997).  Another  important  property  of  the  algorithm  is  that  the 
relative  error  is  actively  controlled  by  prescribing  a  resolution  threshold  parameter.  All  of  these  features 
make  this  algorithm  a  very  attractive  candidate  for  tough  combustion  problems. 
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6  Results 

We  describe  here  two  representative  problems  that  highlight  the  strengths  of  the  WAMR/ILDM  algo¬ 
rithm.  In  the  first  subsection  we  discuss  in  detail  a  one-dimensional  reactive  flow  with  realistic  chemistry. 
In  the  second,  we  give  a  briefer  discussion  of  results  from  a  benchmark  two-dimensional  problem  at  large 
Reynolds  numbers. 

6.1  Ignition  delay  in  a  shock  tube 

We  describe  the  classical  hydrogen-oxygen- argon  shock  tube  ignition  delay  problem.  This  is  a  common 
benchmark  problem  for  both  experiments  and  computations.  Here  we  consider  an  extension  to  the  well- 
documented  case  considered  by  Fedkiw,  et  aL^  1997,  and  fully  described  by  Singh,  et  al.^  2001,  We  consider 
a  shock  tube  of  length  0.12  m  filled  initially  with  H2,  O2,  and  Ar  in  a  2/1/7  molar  ratio.  For  0  <  x  <  0.06  m, 
the  gas  is  taken  to  be  at  p  =  0,18075  kg  m-^,  u  =  487.34  m  andp  =  35.594  kPa,  For  0.06  <  a:  <  0.12  m, 
the  gas  is  at  =  0,072  kg  m”^,  u  =  0  m  s"\  and  p  =  7.173  kPa.  This  state  is  consistent  with  Rankine- 
Hugoniot  jump  conditions  for  the  inviscid  equivalent  of  Equation  (6).  Knowledge  of  these  parameters  allows 
determination  of  all  other  dependent  variables  at  t  =  0  s  through  the  use  of  the  governing  equations.  At 
X  —  0.12  m,  we  consider  a  boundary  which  is  closed  and  adiabatic.  Consequently  u  =  0  m  s  ^  and 
additionally  diffusive  mass  and  energy  fluxes  (n  =  1, . . . ,  A?^),  and  must  be  zero.  At  x  =  0  m,  we 
allow  inflow  conditions  of  u  =  487.34  m  p  —  0.18075  kg  m  ^5  p  =  35.594  kPa,  The  lead  shock  wave  is 
of  insufficient  strength  to  ignite  the  material  on  time  scales  of  interest.  When  the  lead  shock  encounters  a 
wall,  a  stronger  shock  reflects,  and  leaves  the  mixture  at  a  sufficiently  elevated  temperature  so  that  reaction 
soon  ensues.  After  a  brief  ignition  delay,  significant  reaction  commences  at  the  reflecting  wall,  generating 
a  detonation  wave  which  overtakes  the  reflected  shock.  Fedkiw,  et  al.,  1997,  model  this  process  with  one¬ 
dimensional  reactive  Euler  equations,  assuming  a  mixture  of  ideal  gases  with  variable  specific  heats  are 
reacting  in  a  thirty-seven  step  reaction  mechanism  in  the  presence  of  nine  species. 

We  have  first  reproduced  Fedkiw,  et  a/.’s,  1997,  results  using  an  identical  model  and  technique.  It  is 
a  simple  step  to  extend  these  results  to  incorporate  the  ILDM  technique.  We  have  generated  ILDMs  for 
the  entire  range  of  energies  and  densities  encountered  in  this  problem.  Additionally,  we  have  extended  the 
model  to  include  the  effects  of  mass,  momentum  and  energy  diffusion  and  generated  results  using  the  WAMR 
technique  for  spatial  discretization.  Use  of  the  ILDM  closely  replicates  the  full  solution. 

Results  for  the  shock  tube  calculations  are  given  here.  Figure  5  gives  predictions  of  temperature,  velocity, 
pressure,  and  density  vs.  distance  at  time  t  =  195  ps.  At  this  time,  the  lead  right-traveling  inert  shock  has 
reflected  off  the  right  wall  and  is  propagating  to  the  left  with  its  head  near  x  =  0.065  m  at  a  pressure  and 
temperature  of  118.8  kPa  and  1196  K,  respectively.  Close  behind  the  lead  shock  is  the  much  stronger, 
left-propagating  ZND  detonation  wave,  with  its  head  near  x  =  0.072  m.  All  of  the  usual  salient  features  of  a 
ZND  detonation  are  predicted  here.  The  von  Neumann  spike  is  predicted  at  a  pressure  of  around  450  A: Pa, 
and  the  pressure  relaxes  to  near  200  kPa  at  the  right  boundary.  The  post-detonation  temperatures  are  near 
2500  A,  and  the  velocity  is  seen  to  relax  to  a  value  of  zero  at  the  right  boundary. 

The  solid  lines  show  the  predictions  of  the  full  chemical  kinetics  model.  The  dots  show  the  results 
of  the  calculations  using  the  ILDM  resolving  one  reaction  time  scale;  this  can  be  interpreted  as  one-step 
chemistry  with  a  rational  fidelity  to  full  chemical  kinetics.  It  is  seen  on  this  scale  that  the  predictions  are 
nearly  identical!  Examination  of  the  local  eigenvalues  indicates  that  use  of  the  manifold  suppresses  temporal 
resolution  of  reaction  events  which  occur  faster  than  a  roughly  0.1  ps  time  scale.  For  a  given  p  and  e,  we 
construct  the  ILDM  using  adaptive  parameterization  as  described  by  Maas,  1998.  This  is  done  for  sixteen 
values  of  p  ranging  from  0.25  kg  to  1.00  kg  in  steps  of  0.05  kg  m“^.  Similarly  we  use  nineteen  values 
of  e  ranging  from  0.5  x  10^  J  kg~^  to  9.5  x  10^  J  kg~^  in  steps  of  0,5  x  10^  J  kg  ^ .  Hence  304  slices  such  as 
shown  in  Figure  1  were  constructed.  Finally  each  ILDM  was  stored  with  an  equally  spaced  parameterization 
of  100  values  of  Vz/aO  for  easy  lookup.  Thus,  the  ILDM  lookup  table  has  a  size  ofl6xl9xl00. 

The  effects  of  diffusion  are  clearly  seen  when  one  examines  finer  length  scales.  Figure  6  shows  two  views 
of  pressure  vs.  distance  at  a  somewhat  later  time,  t  =  230  ps^  by  which  time  the  detonation  wave  has 
overtaken  the  reflected  shock.  In  the  view  on  the  left,  the  same  length  scale  is  shown  as  in  Fig,  5.  The  view 
on  the  right  shows  a  120  factor  spatial  magnification  near  the  lead  shock.  In  this  figure  the  dots  represent 
the  actual  collocation  points  as  chosen  by  the  WAMR  technique.  It  has  been  found  that  fifteen  levels  of 
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Figure  5:  Predictions  of  temperature,  velocity,  pressure,  and  density  vs.  distance  at  t  =  195  fis  using  a 
maximum  of  300  collocation  points,  15  scale  levels  for  full  chemical  kinetics  (solid  lines)  and  ILDM  kinetics 
(dots)  for  viscous  hydrogen/oxygen/argon  detonation. 


Global  View  Fine  Scale  Structure 


Figure  6:  Predictions  of  pressure  vs.  distance  for  ignition  delay  at  t  =  230  [is  on  coarse  and  fine 

length  scales  demonstrating  the  spatial  resolution  of  viscous  and  induction  zone  structures. 
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wavelets  (equivalent  to  2^^  uniform  grids)  are  sufficient  to  achieve  the  desired  accuracy.  Moreover,  because 
of  the  adaptive  nature  of  the  method,  only  three  hundred  collocation  points  at  most  are  necessary  to  achieve 
this  detail!  It  is  clear  on  this  scale  that  both  the  viscous  shock  and  chemical  induction  zones  have  been 
resolved.  Here  it  is  predicted  that  the  shock  is  essentially  inert  and  has  a  thickness  of  roughly  50  jira. 

The  induction  zone,  a  region  of  essentially  constant  pressure,  temperature,  and  density,  has  a  thickness 
of  roughly  470  fim.  In  the  induction  zone  many  reactions  are  actively  occurring,  giving  rise  to  a  release  of 
energy  which,  because  of  the  extreme  temperature  sensitivity  of  reaction  rates,  accumulates  to  an  extent 
that  a  thermal  explosion  occurs  at  the  end  of  the  induction  zone.  While  the  wavelet  representation  certainly 
has  captured  these  thin  layers,  it  is  noted  that  because  it  was  chosen  not  to  use  individual  species  mass 
fractions  as  part  of  the  adaption  criteria,  some  finer  scale  reaction  zone  structures  have  not  been  spatially 
resolved. 

In  the  process  of  understanding  the  time  scales  associated  with  the  kinetics  of  a  spatially  homogeneous 
mixture,  we  have  computed  all  time  scales  through  an  eigenvalue  analysis.  This  analysis  indicates  that 
reaction  time  scales  as  small  as  sub-nanosecond  are  predicted  by  the  standard  models  of  Maas  and  Warnatz, 
1988,  and  Maas  and  Pope,  1992.  Such  small  time  scales  give  rise  to  small  reaction-induced  spatial  scales 
which  violate  the  continuum  assumption.  It  is  essentially  for  this  reason  that  we  are  not  adapting  our  spatial 
grid  to  capture  the  subsequent  extremely  fine  length  scales  associated  with  individual  species  mass  fraction 
variation.  This  issue  is  pervasive  in  most  calculations  involving  detailed  chemical  kinetics,  but  it  is  not  often 
addressed  since  standard  spatial  discretization  algorithms  are  unable  to  resolve  this  range  of  scales. 

For  this  particular  problem,  use  of  the  full  integration  technique  requires  three  times  as  much  compu¬ 
tational  time  as  the  ILDM  technique.  We  note  however  that  general  conclusions  regarding  computational 
efficiency  are  difficult  to  draw  as  the  savings  realized  will  be  model-dependent  as  well  as  initial  condition- 
dependent.  The  bulk  of  our  savings  are  realized  as  more  and  more  cells  become  chemically  activated.  At  the 
beginning  of  the  calculation,  when  most  cells  are  in  a  cold  state  far  from  equilibrium,  there  is  no  savings. 
The  calculation  itself  took  roughly  ten  hours  on  a  330  MHz  Sun  UltralO  workstation. 

Figure  7  shows  results  for  the  species  mass  fractions  at  t  =  195  Steep  gradients  in  mass  fractions 
are  predicted  near  the  detonation  front.  As  expected,  H02^  H^  and  H2O2  mass  fractions  have  relatively 
small  values  which  peak  at  the  detonation  front.  Under  these  conditions,  the  major  product  is  H2O,  On  the 
length  scales  shown  in  Figs,  5  and  7,  the  results  appear  very  similar  to  the  inviscid  predictions  of  Fedkiw,  et 
a/.,  1997. 

Figure  8  demonstrates  the  adaptive  nature  of  the  WAMR  technique.  It  shows  the  distribution  of  colloca¬ 
tion  points  and  their  scale  levels,  2“-^  ,  j  =  0, . . . ,  14,  at  two  different  times,  first  at  t  =  180  //s,  when  the  lead 
shock  and  the  approaching  detonation  are  present,  and  later  at  t  =  230  /xs,  after  they  have  merged.  In  both 
cases,  at  most  three  hundred  collocation  points  and  fifteen  wavelet  scale  levels  were  sufficient  to  capture  the 
flow  features  to  the  prescribed  dimensionless  error  tolerance  of  10“^.  Moreover,  it  is  clear  that  the  algorithm 
adapts  to  the  features  of  the  flow. 

As  discussed  by  Menikoff,  1994,  inviscid  codes  introduce  pseudo-entropy  layers  near  regions  of  wave- wave 
and  wave-boundary  interactions.  These  often  appear  as  0(1)  anomalies  in  temperature  and  density  near  the 
wall  in  shock  tube  predictions.  Figure  9  shows  the  results  of  our  viscous  calculation  in  a  spatial  zone  near 
the  wall  just  after  shock  reflection.  On  this  scale,  there  is  no  apparent  entropy  layer  near  the  wall. 

A  finer  scale  examination  of  the  dependent  variables,  shown  in  Figure  10,  reveals  what  is  happening. 
It  is  evident  from  the  temperature  plot  that  there  is  a  small  entropy  layer  near  the  wall,  here  physically 
induced.  The  physical  diffusion  mechanisms  rapidly  smear  the  layer  within  a  few  microseconds.  It  may  be 
possible  that  the  correct  capturing  of  a  temperature-sensitive  ignition  event  near  a  wall  could  be  critically 
dependent  on  having  the  correct  physics  in  the  model.  For  the  viscous  calculation,  a  temperature  rise  of 
roughly  6  K  is  predicted.  Performing  a  similar  calculation  with  an  inviscid  Godunov-based  model  with 
first  order  upwind  spatial  discretization  near  regions  of  steep  gradients  using  400  evenly  spaced  grid  points 
induces  a  temperature  rise  of  nearly  20  K ,  which  persists.  It  might  be  expected  that  numerical  diffusion 
would  dissipate  this  temperature  spike.  However,  as  detailed  by  LeVeque,  1992,  the  leading  order  numerical 
diffusion  coefficient  for  such  methods  is  proportional  to  the  local  fluid  particle  velocity.  As  the  fluid  particle 
velocity  at  the  wall  and  in  the  region  downstream  of  the  reflected  shock  is  zero,  the  effects  of  numerical 
diffusion  here  are  at  most,  confined  to  higher  order.  This  temperature  rise,  similar  to  that  obtained  by 
Fedkiw  et  al,  1997,  1999,  has  been  obtained  by  effectively  imposing  an  adiabatic  boundary  condition  done 
through  extrapolation.  This  type  of  condition  is  inconsistent  with  the  inviscid  governing  equations. 
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Figure  7:  Predictions  of  species  mass  fractions  vs.  distance  at  t  =  195  fis  using  a  maximum  of  300  col¬ 
location  points,  15  scale  levels  for  full  chemical  kinetics  (solid  lines)  and  ILDM  kinetics  (dots)  for  viscous 
hydrogen /oxygen /argon  detonation. 


Figure  8:  Spatial  distribution  of  collocation  points  and  levels  at  ^  =  180  (two  shock  structure)  t  —  230  fis 
(single  shock  structure)  demonstrating  grid  adaption. 
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Figure  9:  Predictions  of  temperature,  velocity,  pressure,  and  density  vs.  distance  before  commencement  of 
significant  reaction  but  after  shock  reflection  using  a  reactive  Navier-Stokes  model,  t  ~  177  fis. 


X  (m)  X  (m) 

Figure  10:  Close-up  view  of  predictions  of  temperature,  velocity,  pressure,  and  density  vs.  distance  before 
commencement  of  significant  reaction  but  just  after  shock  reflection  {t  =  174  fis)  and  slightly  later  {t  = 
177  fis). 
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6.2  Lid-driven  cavity 

Here  results  are  presented  for  a  standard  benchmark  flow  using  the  WAMR  algorithm:  a  lid-driven  fluid  in  a 
square  cavity.  In  this  problem,  a  fluid  is  confined  in  a  square  cavity.  Three  of  the  walls  are  stationary,  while 
one  wall  has  a  constant  velocity.  In  order  to  satisfy  the  no-slip  boundary  condition  at  all  surfaces,  the  fluid 
is  driven  into  rotational  motion.  As  the  Reynolds  number  is  increased,  more  and  more  complex  structures 
of  different  scales  are  generated  in  the  flow  field;  such  a  flow  is  thus  a  good  candidate  for  study  with  the 
WAMR  technique.  We  compare  our  results  with  the  widely  cited  study  of  Ghia,  et  al.^  1982,  who  considered 
a  finite  difference  solution  of  precisely  the  same  problem  using  multigrid  methods  in  which  the  finest  grid 
was  equivalent  to  a  257  x  257  uniformly  spaced  computational  grid. 

The  dimensionless,  unsteady,  two-dimensional,  incompressible  Navier-Stokes  equations  are  adopted  as  a 
model: 


du 

'm 

dv 


du  dv 

dx  dy 

du  du 

dx  dy 

dv  dv 
+  -bu— 
dx  dy 


dp  1  /  d'^u  d'^u\ 

dx  Re  ^  dy"^ )  ' 
dp  \  ( d'^v 

dy  ^  Re  ^  dy'^  )  ’ 


(28) 

(29) 

(30) 


In  these  equations  independent  variables  are  Cartesian  coordinates  x  and  y,  and  time  t.  Dependent  variables 
are  rc-velocity  u,  y-velocity  u,  and  pressure  p.  The  only  parameter  is  the  Reynolds  number,  Re, 

In  applying  the  WAMR  technique  to  this  problem,  the  dependent  variables  are  projected  onto  a  wavelet 
basis  extended  to  two-dimensions.  Many  of  the  same  techniques  discussed  for  one-dimensional  methods  apply 
here  as  well,  and  will  not  be  discussed.  Because  of  the  two-dimensional  nature  of  the  flow  however,  it  was 
necessary  to  develop  a  robust  computational  algorithm  for  computing  derivatives  on  a  non-uniform  stencil. 
Further,  a  fast  Poisson  solver  was  developed  to  solve  for  the  elliptic  equation  for  pressure  that  was  necessary 
to  solve  at  every  time  step.  Dynamic  memory  allocation  was  utilized  for  fluctuating  storage  requirements. 

Results  at  Re  =  10000,  the  highest  Reynolds  number  in  the  published  literature  for  which  results  exist, 
are  shown  in  Figures  11,  12,  and  13,  at  early  stages  of  the  flow  development,  at  t  =  4, 12, 100,  respectively. 
It  is  noted  that  the  asymptotic  flow  is  oscillatory.  Each  of  these  figures  gives  a  plot  of  the  velocity  vector 
field.  It  is  important  here  to  highlight  the  transitory  portion  of  the  development  so  that  the  adaptive  nature 
of  the  method  can  be  observed.  In  these  figures,  a  velocity  vector  is  affixed  to  the  location  of  each  wavelet. 
Figure  11  shows  a  prominent  starting  vortex,  and  there  is  a  high  concentration  of  wavelets  in  regions  of  steep 
gradients  of  the  flow.  In  all  of  the  calculations  reported,  five  wavelet  levels  have  been  utilized.  The  coarsest 
set  of  wavelets  gives  results  equivalent  to  a  16  x  16  grid  and  the  finest  is  equivalent  to  at  256  x  256  grid,  and 
is  an  equivalent  resolution  as  Ghia,  et  a/.’s  finest  grid.  Figure  12  gives  a  snapshot  at  a  later  time  when  the 
primary  vortex  has  convected  and  diffused  further  into  the  cavity.  Once  again,  it  is  obvious  that  the  wavelet 
collocation  points  have  adapted  to  the  changing  nature  of  the  flow.  Figure  13  shows  results  at  a  much  later 
time.  Here,  careful  examination  (not  shown  here)  reveals  that  recirculation  zones  in  the  corners  have  been 
predicted. 

Results  at  a  lower  Reynolds  number,  Re  =  3200,  are  shown  in  Figure  14.  Because  of  the  lower  Reynolds 
number,  the  fluid  has  relaxed  to  a  steady  state  at  the  depicted  time  of  i  =  100.  Figure  15  depicts  the 
u  component  of  velocity  as  a  function  of  y  for  2:  =  0.5.  The  results  from  the  WAMR  are  in  very  good 
agreement  with  published  results  of  Ghia,  et  ai  Moreover,  we  believe  our  results  are  actually  more  accurate 
since  Ghia  et  ai  found  it  necessary  to  use  upwind  differencing  for  convective  terms  to  stabilize  the  numerical 
computations.  As  is  well  known,  the  use  of  upwind  differencing  introduces  numerical  diffusion  into  the 
solution.  Figure  16,  similar  to  Figure  15,  depicts  the  v  component  of  velocity  as  a  function  of  2:  for  y  =  0.5. 
Again  the  results  are  very  close  to  those  of  Ghia,  et  ai  at  the  scale  shown  in  the  figure. 


7  Conclusion 

The  essential  goals  of  this  research  project  have  been  achieved.  Significant  progress  has  been  made  in 
the  development  of  both  the  method  of  intrinsic  low  dimensional  manifolds  as  well  as  the  wavelet  adaptive 
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Ra  =  10000,  time  a  4 


Figure  11:  Velocity  vector  field  and  collocation  points  for  lid-driven  cavity,  Re  —  10000,  t  —  4. 


Re  =10000,  lime  =12 


Figure  12:  Velocity  vector  field  and  collocation  points  for  lid-driven  cavity,  Re  —  10000,  t  —  12. 


multilevel  representation  algorithms.  The  methods  have  successfully  been  employed  on  challenging  multiscale 
problems  in  reactive  and  multidimensional  inert  flows  to  predict  results  which  have  met  or  exceeded  existing 
published  standards. 

The  work  is  not  incremental  in  character  and  in  fact  represents  a  substantially  new  approach  which  we 
believe  has  yielded  consequent  dividends.  While  additional  generalizations  can  and  should  be  made  to  the 
algorithm,  we  nevertheless  believe  the  research  has  set  a  new  standard  in  high  performance  computing  applied 
to  combustion.  In  addition  we  believe  the  combined  novelty  and  utility  of  the  WAMR/ILDM  method  will 
have  ramifications  extending  into  the  other  areas  of  computational  physics. 


20 


Re=  10000,  lime=  100 


Figure  13:  Velocity  vector  field  and  collocation  points  for  lid-driven  cavity,  Re 


10000,  t 


Re  =  3200,  time  =  100 


=  100. 


Figure  14:  Velocity  vector  field  and  collocation  points  for  lid-driven  cavity,  Re  =  3200,  steady  state  results 
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Figure  15:  Comparison  of  midplane  (x  =  1/2)  u  with  Ghia,  et  aZ.’s,  1982,  results  for  lid-driven  cavity  at 
Re  =  3200. 


V-vdodty  along  y  =  0.5,  Re  =  3200 


x-coordinate 


Figure  16:  Comparison  of  midplane  {y  =  1/2)  v  with  Ghia,  et  a/.’s,  1982,  results  for  lid-driven  cavity  at 
Re  =  3200 
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